High-density amorphous ice: A path-integral simulation 
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Structural and thermodynamic properties of high-density amorphous (HDA) ice have been studied 
by path-integral molecular dynamics simulations in the isothermal-isobaric ensemble. Interatomic 
interactions were modeled by using the effective q-TIP4P/F potential for flexible water. Quantum 
nuclear motion is found to affect several observable properties of the amorphous solid. At low tem- 
perature (T = 50 K) the molar volume of HDA ice is found to increase by 6%, and the intramolecular 
O-H distance rises by 1.4% due to quantum motion. Peaks in the radial distribution function of 
HDA ice are broadened respect to their classical expectancy. The bulk modulus, B, is found to 
rise linearly with the pressure, with a slope dB/dP = 7.1. Our results are compared with those 
derived earlier from classical and path-integral simulations of HDA ice. We discuss similarities and 
discrepancies with those earlier simulations. 

PACS numbers: 61.43.Er, 65.60.+a, 62.50.-p, 71.15.Pd 
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I. INTRODUCTION 

Water is ubiquitous on Earth, and a detailed under- 
standing of its properties and behavior under different 
conditions is of crucial importance in several scientific 
fields^— This does not only refer to thermodynami- 
cally stable water phases, but also to various metastablc 
phases obtained at ambient and extreme conditions of 
temperature and pressure. In particular, several forms 
of amorphous ice have been found and studied by both 
experimental^— and theoretical^— methods, but some 
of their properties still lack a complete understanding. 
This is mainly due to the peculiar structure of condensed 
phases of water, where hydrogen bonds between adjacent 
molecules give rise to rather unusual properties of these 
phases. 

The atomic dynamics in amorphous solids cause the 
appearance of localized low-energy excitations, display- 
ing appreciable deviations from the situation of atomic 
nuclei harmonically vibrating around their potential 
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Such deviations from harmonicity, com- 
bined with the quantum character of atomic dynamics, 
is of great importance in the characterization of amor- 
phous materials. In this context, the classical papers by 
Phillips^ 7 - and Anderson et al*£ opened a research line 
by modeling low-energy excitations in amorphous solids 
by two-level systems. In more recent years there have 
appeared several detailed descriptions of the low-energy 
motion in this type of materials, beyond the standard 
tunneling models 

When one studies average structural properties of 
amorphous materials, an interesting question is whether 
quantum effects can be noticeable in the presence of 
structural disorder. To be concrete, one may ask whether 
the radial distribution function (RDF) obtained in clas- 
sical simulations is appreciably modified by considering 
quantum atomic derealization, which can be important 
mainly at low temperatures. Two factors appear to com- 
pete in the broadening of the RDF peaks at low temper- 



ature: structural disorder and quantum derealization. 20 
As a first approach, one expects that for amorphous ma- 
terials with heavy atoms (small zero-point vibrational 
amplitudes), structural disorder will broaden the peaks 
more than zero-point motion, and the opposite may occur 
for disordered materials with light atoms. For amorphous 
ice, one may suspect that the presence of hydrogen will 
make quantum effects appreciable, even for strong struc- 
tural disorder. 

Computer modeling of amorphous ice has been em- 
ployed in recent years to obtain insight into its struc- 
tural and dynamical properties . 12 i 13 ] 21 ~ — The beginning 
of computer simulations of condensed phases of water 
at an atomic level dates back more than 40 yearSj 25 i 26 
and nowadays a large variety of empirical interatomic 
potentials can be found in the literature.— ~— Many of 
them assume a rigid geometry for the water molecule, 
whereas some others allow for molecular flexibility ei- 
ther with harmonic or anharmonic OH stretches. In re- 
cent years, simulations of water using ab initio density 
functional theory (DFT) have been also carried out.— ~— 
However, hydrogen bonds in condensed phases of water 
are difficult to describe with currently available energy 
functionals, making that some properties are not accu- 
rately reproduced by DFT calculations, 37 Some progress 
in the description of van der Waals interactions in water 
within the DFT formalism has been recently made<22r— 

A shortcoming of ab-initio electronic-structure calcu- 
lations is that they usually deal with atomic nuclei as 
classical particles, disregarding quantum effects like zero- 
point motion. These effects may be accounted for using 
harmonic or quasiharmonic approximations for the nu- 
clear motion, but the precision of these approaches is not 
readily estimated when large anharmonicities are present, 
as can be the case for light atoms like hydrogen in dis- 
ordered materials. To take into account the quantum 
character of atomic nuclei, the path-integral molecular 
dynamics (PIMD) approach has proved to be very useful, 
since in this procedure the nuclear degrees of freedom can 
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be quantized in an efficient manner, thus including both 
quantum and thermal fluctuations in many-body systems 
at finite temperatures.— This computational technique is 
now well established as a tool to study problems in which 
anharmonic effects can be important^ 3 - Thus, a powerful 
approach could consist in combining DFT to determine 
the electronic structure and path integrals to describe the 
quantum motion of atomic nuclei However, this pro- 
cedure requires computer resources that would restrict 
enormously the number of state points that can be con- 
sidered in actual calculations. 

Several forms of amorphous ice have been de- 
tected in recent years, corresponding to different 
densitieSi 12 i 14 ' 44 " — In the present paper we study high- 
density amorphous (HDA) ice by PIMD simulations at 
different pressures and temperatures, to analyze its struc- 
tural and thermodynamic properties. Interatomic inter- 
actions are described by the flexible q-TIP4P/F model, 
which was recently developed and has been employed to 
carry out PIMD simulations of liquic— — and solid^Sr— 
water. Here we pose the question of how quantum motion 
of the lightest atom can influence the structural proper- 
ties of an amorphous water phase, and in particular if this 
quantum motion is appreciable for the solid at different 
densities, i.e. under different external pressures. This 
refers to the crystal volume and interatomic distances, 
but also to the mechanical stability of the solid. These 
questions have been addressed earlier for high- and low- 
density amorphous icei 21 ' 53 In this context, it is usually 
assumed that increasing quantum fluctuations enhances 
the exploration of the configuration space, but in certain 
regimes an increase in quantum fluctuations can lead to 
dynamical arrest, as found for glass formation^ 

The paper is organized as follows. In Sec. II, we de- 
scribe the computational method and the model em- 
ployed in our calculations. In Sec. Ill we discuss the 
method employed to generate the simulation cells of 
amorphous ice. Our results are presented in Sec. IV, 
dealing with molar volume, interatomic distances, ki- 
netic energy, and bulk modulus of HDA ice. Sec. V gives 
a comparison with earlier work, and Sec. VI includes a 
summary of the main results. 



II. COMPUTATIONAL METHOD 

We employ PIMD simulations to obtain several proper- 
ties of amorphous ice at different temperatures and pres- 
sures. This kind of simulations are based on an isomor- 
phism between a quantum system and a classical one, 
obtained after a discretization of the quantum density 
matrix along cyclic paths . 55 ' 56 This isomorphism corre- 
sponds to replacing each quantum particle by a ring poly- 
mer consisting of L (Trotter number) classical particles, 
joined by harmonic springs with temperature- and mass- 
dependent force constant. Details on this simulation pro- 
cedure can be found elsewhere . 42 ' 43 The dynamics used 
in this method is fictitious and does not correspond to 



the real quantum dynamics of the considered particles, 
but it helps to effectively sample the many-body config- 
uration space, yielding precise results for the properties 
of the actual quantum system. Another way to derive 
such properties could be the use of Monte Carlo sam- 
pling, but we have found that this procedure requires 
for the present problem more computer resources than 
PIMD simulations. An important advantage of the lat- 
ter is that the computing codes can be more readily par- 
allelized, which turns out to be a relevant factor for an 
efficient use of modern computer architectures. 

Simulations of crystalline and amorphous ice were car- 
ried out here in the isothermal-isobaric NPT ensemble 
(TV, number of particles; P, pressure; T, temperature), 
which allows one to find the equilibrium volume of a solid 
at given pressure and temperature. We have used ef- 
fective algorithms for carrying out PIMD simulations in 
this statistical ensemble, similar to those described in 
the literature; 5 -^— and employed earlier in the study of 
solid and liquid water by PIMD simulations. We have 
considered temperatures between 50 K and 300 K, and 
pressures up to 8 GPa. Both negative (tension) and pos- 
itive (compression) pressures have been employed in the 
simulations. For negative P we considered amorphous ice 
in the region of mechanical stability of the solid, down 
to P ~ — 0.5 GPa. For comparison with results of PIMD 
simulations, some simulations of HDA ice have been also 
performed in the classical limit, which is obtained in our 
path integral procedure by setting the Trotter number 
L = 1. 

Our PIMD simulations were carried out on cells in- 
cluding 96 or 216 water molecules, which were generated 
from pressure treatment of ice Ih and Ic supercells, re- 
spectively. The former (96 molecules) corresponds to a 
(3a, 2y/3a, 2c) supercell, where a and c are the standard 
hexagonal lattice parameters of ice Ih, whereas the lat- 
ter (216 molecules) corresponds to a 3 x 3 x 3 supercell 
of the cubic unit cell of ice Ic. The main purpose of 
taking these two ice types was to check the influence of 
the starting crystalline ice on the properties of the amor- 
phous ice resulting under pressure. For all considered 
variables (structural and thermodynamic), we found in 
both cases results that coincided with each other within 
statistical error bars of the simulation procedure. To gen- 
erate proton-disordered ice supercells (of Ih and Ic types) 
prior to amorphization, we employed a Monte Carlo pro- 
cedure to impose that each oxygen atom had two chem- 
ically bonded and two H-bonded hydrogen atoms, and 
with a cell dipole moment close to zeroi 51 ' 60 For some par- 
ticular conditions, we checked that both starting proton 
disorder and history of the amorphous supercell (amor- 
phization procedure) do not affect significantly the re- 
sults presented below. In particular, once formed the 
amorphous material from crystalline ice, we checked that 
its properties are reversible upon increasing and decreas- 
ing temperature and/or pressure. 

Interatomic interactions were modeled by the point 
charge, flexible q-TIP4P/F model, developed to study 



liquid water ^ and that was later employed to study var- 
ious properties of ic o 50 ' 51 and water clusters d Many of 
the empirical potentials used earlier for quantum simula- 
tions of condensed phases of water treat H2O molecules 
as rigid bodies^r— This turns out to be convenient 
for computational efficiency, but neglects the role of in- 
tramolecular flexibility in the structure, dynamics, and 
thermodynamics of the condensed water phases4£> More- 
over, the q-TIP4P/F potential takes into account the sig- 
nificant anharmonicity of the O-H vibration in a water 
molecule by considering anharmonic stretches, vs. the 
harmonic potentials employed in most of the simulations 
that considered quantum effects in these water phases. 

Technical details on the PIMD simulations presented 
here are the same as those described and used in 
Refs. [50ll5ll The Trotter number L has been taken pro- 
portional to the inverse temperature (L oc 1/T), so that 
LT = 6000 K, which allows us to keep roughly a constant 
precision in the PIMD results at the different tempera- 
tures under consideration. The time step At employed 
for the calculation of interatomic forces in the molecular 
dynamics procedure was taken in the range between 0.1 
and 0.3 fs, which was found to give adequate convergence 
for the variables studied here. For given conditions of 
pressure and temperature, a typical simulation run con- 
sisted of 2 x 10 5 PIMD steps for system equilibration, 
followed by 10 6 steps for the calculation of ensemble av- 
erage properties. 
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FIG. 1: Molar volume of ice as a function of pressure, as 
derived from PIMD simulations at T — 75 K. Open circles 
represent results derived from simulations starting from ice 
Ih. Solid circles are data points obtained from simulations 
starting from amorphous ice. Error bars are in the order 
of the symbol size. The dashed line is a guide to the eye. 
An open triangle at P = indicates the volume measured 
by Mishima et al^ The solid line was obtained from the 
pressure-density data displayed in the review by Loerting and 
Giovambattistai^, taken from Ref. Ir36l . 



III. PREPARATION OF AMORPHOUS ICE 

In this Section we give details on the form in which we 
obtain the disordered structures of HDA ice that are sub- 
sequently employed to characterize and study this amor- 
phous phase from PIMD simulations. In particular, we 
have obtained simulation cells of HDA ice by applying a 
hydrostatic pressure to ice Ih and ice Ic at several tem- 
peratures. This procedure allows us also to check the 
pressure at which amorphization occurs, and to compare 
it with data (both experimental and theoretical) reported 
in the literature. We note that HDA ice has been recently 
obtained at room temperature from ice VII under rapid 
compression^ 

In Fig. 1 we present the pressure dependence of the 
molar volume of ice at a temperature of 75 K, as derived 
from our PIMD simulations. Results shown as open cir- 
cles correspond to simulations starting from ice Ih, for 
both negative and positive pressures. We observe that in 
the region between -1 GPa and 1 GPa the volume de- 
creases smoothly as pressure is increased, and at about 
1.2 GPa it suffers a sudden decrease which corresponds 
to ice amorphization. This value is close to the spinodal 
pressure (limit of mechanical stability) obtained for ice 
Ih at this temperature in Ref. 67 (P s = 1.19±0.05 GPah 
and to the amorphization pressure obtained in Ref. [22j 
from classical molecular dynamics simulations. From 1 
to 2 GPa the volume decreases by 27%, and at higher 



pressures it continues decreasing to reach a value of 10.8 
cm 3 /mol at P = 8 GPa, to be compared with v = 19.36 
cm 3 /mol found for ice Ih at atmospheric pressure. 

We have carried out this procedure at several tem- 
peratures, and the results are similar to those shown in 
Fig. 1. In particular, the amorphization pressure changes 
slightly with temperature, and at 250 K we find 0.95 
GPa. This was discussed in Ref. [67| in connection with 
the stability of ice Ih under pressure, and will not be 
repeated here. Apart from the pressure at which the 
solid amorphizes in the simulations, more precise values 
for the amorphization pressure can be obtained from the 
spinodal line of ice Ih, which gives the limit for the me- 
chanical stability of this solid phase. This pressure corre- 
sponds to the vanishing of the bulk modulus (divergence 
of the compressibility) , and can be approached in com- 
puter simulations at low temperatures (T < 50 K). At 
higher T ice Ih amorphizes in the PIMD simulations be- 
fore reaching the corresponding spinodal pressure, due to 
nucleation effects leading to the breakdown of the ice Ih 
structured We note that the formation of HDA ice from 
ice Ih has not been observed in the laboratory at tem- 
peratures so high as 250 K, but we obtain this transition 
at such temperatures and find an amorphous phase that 
remains metastable along our simulations (which in fact 
can only cover time scales shorter than those of actual ex- 
periments). HDA ice has been, however, obtained from 
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FIG. 2: Oxygen-oxygen radial distribution function of ice Ih 
(dashed line) and HDA ice (solid line) at T = 100 K and P 
= 1 atm, as derived from PIMD simulations. 



solid remains close to that of ice Ih, in spite of the appre- 
ciable volume reduction suffered by the solid in the anror- 
phization process. Note that the difference in height of 
the first peaks in the crystalline and amorphous phases is 
mainly due to the definition of radial distribution func- 
tion, which is normalized by the density^ In fact, we 
calculate g(r) as 



50) = 



d(r) 
D 



(1) 



where D is the mean density (atoms per unit volume) and 
d(r) is the "local density" at distance r from a reference 
atom: 



d(r) 



N(r) 
4-7rr 2 Ar ' 



(2) 



N(r) being the number of atoms at distances between r 
and r + Ar. 



IV. PROPERTIES OF HIGH-DENSITY 
AMORPHOUS ICE 



ice VII at room temperature^ as mentioned above. 

Solid water amorphizes in an irreversible way, so that 
new PIMD simulations at pressures lower than 1 GPa, 
starting from the amorphous phase, do not recover the 
crystalline phase. This is shown in Fig. 1 as solid sym- 
bols. In these simulations the pressure was gradually re- 
duced down to negative pressures until reaching the limit 
of mechanical stability of the material. At 75 K we could 
reach a pressure of -0.4 GPa, and at still more negative 
pressures the amorphous solid broke down along the sim- 
ulations, transforming into the gas phase. This point will 
be discussed below in connection with the bulk modulus 
of HDA ice. For comparison with our results we present 
in Fig. 1 the molar volume obtained by Mishima et al— 
for HDA ice from x-ray diffraction experiments at 80 K 
and atmospheric pressure. Shown is also the volume- 
pressure curve derived from the data given in the review 
by Loerting and Giovambattista^ taken from Ref. [66| 
(solid line). We note that between P = and 1 GPa 
the molar volume of the amorphous material is clearly 
smaller than that of the crystal, but at negative pressure 
the volume of the amorph increases fast, due to the prox- 
imity of its metastability limit, and therefore approaches 
the volume of ice Ih. 

For additional confirmation of the amorphous charac- 
ter of the solid obtained after the cycle indicated by ar- 
rows in Fig. 1 (first pressure increase up to 8 GPa, 
and then pressure reduction), we display in Fig. 2 the 
0-0 RDF for ice Ih and amorphous ice, as derived from 
PIMD simulations at T = 100 K and atmospheric pres- 
sure. Upon amorphization, the prominent peaks corre- 
sponding to the second and third coordination sphere 
merge into one broad feature with a maximum at about 
4.1 A. The maximum of the first peak in the amorphous 



A. Volume 

As shown above in Fig. 1 and discussed in Sect. Ill, ice 
Ih suffers an important reduction in volume upon amor- 
phization at about 1.2 GPa. This volume decrease is 
associated to a softening of intermolecular O-H bridges, 
accompanied by a reduction in the mean distance to oxy- 
gen atoms in the second and third coordination shells (see 
Fig. 2). After amorphization, the molar volume contin- 
ues decreasing as pressure is raised, and we find at 75 
K a reduction from 13.11 cm 3 /mol at P = 2 GPa to 
10.81 cm 3 /mol at 8 GPa, which means a decrease of 18% 
respect to the volume at 2 GPa. 

We emphasize that the volume of the amorph decreases 
smoothly in the whole pressure region considered here, 
and we did not detect in this region any other phase 
change as those reported in the literature. In fact, Hem- 
ley et alX observed a pressure-induced re-crystallization 
of HDA ice at about 4 GPa. Also, a slow transformation 
of HDA to cubic ice on slow depressurization has been 
observed^ In our simulations the amorph remains in its 
metastable state in the whole range of temperature and 
pressure studied here. 

In Fig. 3 we show the molar volume of ice Ih (circles) 
and the HDA phase (squares) as a function of tempera- 
ture at atmospheric pressure, as derived from our PIMD 
simulations. At P = 1 atm and T = 75 K, after re- 
leasing the pressure applied for amorphization, we find a 
molar volume v — 15.57 cm 3 /mol, which corresponds to 
a density p — 1.16 g/cm 3 , in good agreement with the 
values given by Mishima et al. at zero pressure: p — 
1.17 g/cm 3 in Ref. [H and 1.19 g/cm 3 in Ref. 0. For com- 
parison with the PIMD results, we also display in Fig. 3 
classical data derived here with the q-TIP4P /F potential 
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FIG. 3: Temperature dependence of the molar volume of ice 
at atmospheric pressure as derived from PIMD simulations: 
Ih (circles) and HDA ice (squares). Results of classical sim- 
ulations for HDA ice are displayed as triangles. A solid line 
represents data obtained from classical molecular dynamics 
simulations by Tse et al*^ For ice Ih the error bars are smaller 
than the symbol size. Lines are guides to the eye. 



(triangles), as well as those obtained by Tse et alM us- 
ing the SPC/E potential. At low temperature, the molar 
volume found by these authors is similar to that found 
here in the classical simulations, and it becomes closer to 
the PIMD results as temperature rises. Comparing our 
results with the q-TIP4P/F force field, we find at 50 K 
an increase in molar volume of about 0.85 cm 3 /mol due 
to nuclear quantum effects, which amounts to a 6% of 
the classical value. 

In the temperature region up to 250 K we did not ob- 
serve any transition from HDA ice to a low-density amor- 
phous phase with density p = 0.94 g/cm 3 (molar volume: 
19.1 cm 3 /mol), as experimentally observed at about 120- 
130 Kj£&2iii and HDA ice remained as a metastable 
phase along our PIMD simulation runs. Something sim- 
ilar happens with the classical simulations reported in 
Ref. [Tj. We believe that such a transition between amor- 
phous phases is not captured by the simulations due to 
the short time window that in fact can be observed in 
the calculations, similarly to the difficulties found in this 
kind of simulations to directly obtain transitions between 
different crystalline phases. This question could possi- 
bly be solved by performing direct coexistence simula- 
tions as those reported in the literature for liquid-solid 
transitions.— 

It is clear that the amorphous material has a higher 
density, or a smaller molar volume than ice Ih, but it ex- 
pands with temperature much faster than the crystalline 
solid. At atmospheric pressure we find for the amorph a 
volume increase of 2.6 cm 3 /mol in the temperature range 



0.985 



0.98 

U 
O 

d 

B 0.975 

X 

o 



0.965 



0.96 





1 I 


1 1 1 1 1 1 




T = 75 K 














_»-—-—■ 


quantum 






. • • ' ' 










classical 

i , i > i i i 





2 


4 6 8 



Pressure (GPa) 

FIG. 4: Mean intramolecular O-H distance as a function of 
pressure for amorphous ice at 75 K, as derived from PIMD 
(squares) and classical (circles) simulations. Error bars are in 
the order of the symbol size. Lines are guides to the eye. 



from 50 to 250 K. In the same temperature region, the 
volume of ice Ih is found to rise less than 0.3 cm 3 /mol, in 
part due to the negative thermal expansion around 100 
K (not clearly observed at the scale of Fig. 3). Thus, we 
find for the thermal expansion coefficient of HDA ice at 
100 K: a = 8.1 x 10~ 4 K _1 . This large volume expansion 
of the amorph, as compared with ice Ih, is indeed related 
with the lower bulk modulus of the amorphous material 
at P = 1 atm (see below). 



B. Interatomic distances 

In this section we present results for interatomic dis- 
tances in amorphous ice between atoms in the same and 
adjacent molecules. This can shed light on the struc- 
tural changes suffered by the material when temperature 
and/or pressure are modified. We first show in Fig. 4 the 
mean intramolecular O-H distance as a function of pres- 
sure at a temperature T — 75 K. In connection with this 
figure, there are two results that should be emphasized. 
First, at a given pressure, we find that the O-H bond 
distance increases appreciably due to nuclear quantum 
effects. In fact, a classical simulation at T — 75 K and 
ambient pressure yields a mean O-H distance of 0.966 A, 
to be compared with 0.980 A derived from PIMD simula- 
tions, which means an increase of 1.4% in the bond length 
due to nuclear quantum motion. This difference is much 
larger than the temperature-induced change in d(O-H) 
at atmospheric pressure, which amounts to about 0.002 
A in the range from 50 to 250 K. The increase due to 
quantum motion is rather constant in the whole pressure 
range studied here. 



Another important result observed in Fig. 4 is that at 
a given temperature the O-H distance increases as pres- 
sure is raised, contrary to the usual contraction of atomic 
bonds for increasing P. This somewhat anomalous fact 
is due to the characteristic structure of ice with hydro- 
gen bonds connecting contiguous water molecules, which 
gives rise to an anticorrelation between the strength of 
molecular O-H bonds and intermolecular H bridges^ In 
fact, increasing the pressure causes a hardening of in- 
termolecular H bridges, with an associated weakening of 
intramolecular O-H bonds, and a concomitant increase 
in the bond length. This weakening of intramolecular 
O-H bonds in ice for rising pressure has been observed 
experimentally and reported in the literature^ 

It is interesting to compare the O-H bond distance in 
amorphous ice with that obtained for ice Ih in the same 
type of PIMD simulations. At atmospheric pressure and 
75 K, we hnd for ice Ih a mean distance d(O-H) = 0.984 
A, i.e. about 0.4% longer than in HDA ice at the same 
conditions^ At the same temperature and P = 1 GPa, 
close to the amorphization pressure of ice Ih we found a 
distance difference of 0.6%. This means that the average 
O-H distance decreases upon amorphization of the solid, 
which is consistent with a weakening of H bridges be- 
tween adjacent molecules, causing a strengthening of the 
intramolecular bonds, and therefore a shortening of the 
corresponding O-H distance. These differences between 
O-H distances in Ih and HDA ice are consistent with 
changes in the stretching vibrational frequencies, as ob- 
served from Raman scattering experiments. In fact, for 
HDA ice one observes at 80 K a broad Raman band with 
a maximum at rs 3200 cnr 1 ^ to be compared with the 
largest feature appearing in the O-H stretching region 
of ice Ih at about 3090 cm" 1 ^ This hardening of the 
stretching vibrations upon amorphization is consistent 
with the general trend found for water molecules, when 
the intramolecular O-H distance decreases in different 
crystal surroundings^ Aw(0-H)/Ad(0-H) = -2.4 xlO 4 
cm _1 /A. We note, for comparison, that Bellissent-Funel 
et a/.— found an intramolecular O-D distance of 0.97 
A, from neutron scattering experiments on high-density 
amorphous D2O. 

Another related aspect of the O-H distance is its tem- 
perature dependence. For ice Ih at atmospheric pres- 
sure, this distance is known to decrease as temperature is 
raised, as a consequence of the hardening of the bond for 
increasing temperature. This is in line with an enhance- 
ment of molecular motion for rising T, which causes a 
weakening of the H bridges and an associated enhance- 
ment of intramolecular bond strength. Something simi- 
lar is found for amorphous ice in our PIMD simulations, 
where the average O-H distance decreases from 0.9804(2) 
to 0.9788(2) A when T rises from 75 to 250 K. 

Given that nuclear quantum motion is appreciable in 
the mean O-H distance in water molecules in amorphous 
ice, it is expected that isotopic effects can be observed in 
the RDF. In Fig. 5 we display the O-H radial distribution 
functions for HDA ice, as derived from PIMD simulations 
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FIG. 5: Oxygen-hydrogen radial distribution function at 75 
K and P = 1 atm, as derived from quantum PIMD simula- 
tions for H2O (solid line) and D2O (dashed line) amorphous 
ice, as well as from classical molecular dynamics simulations 
(dashed-dotted line). Inset: RDF in the region around 1 
A, showing the peak corresponding to intramolecular O-H 
bonds. 



for normal and deuterated water, as well as from classical 
molecular dynamics simulations. The RDF derived from 
classical simulations is similar to that found earlier in 
this kind of simulations for HDA iceJ^ For interatomic 
distances r from 1.5 to 5 A, we observe in Fig. 5 that 
quantum effects cause a broadening of the peaks. This 
is particularly observable in the peaks at about 1.8 A 
and 3.2 A. The first peak in this RDF, corresponding 
to the intramolecular O-H bonds, is much higher, and 
is displayed in the inset. For the classical model, it is 
about 15 times larger than the peak at 1.8 A. All peaks 
are found to broaden due to quantum effects, and their 
widths are larger for smaller isotopic mass. In fact, for 
the peak corresponding to intramolecular O-H bonds, we 
obtained a full width at half maximum of 0.05, 0.14, and 
0.16 A, for classical ice, quantum D2O, and quantum 
H2O, respectively. Note that in this respect the classical 
limit behaves as the large-mass limit. 

In Fig. 6 we show the hydrogen- hydrogen RDF of HDA 
ice as derived from classical (dashed line) and PIMD sim- 
ulations (solid line) . As expected, quantum motion of H 
broadens the peaks in the RDF with respect to the clas- 
sical result. The peak corresponding to H-H pairs inside 
water molecules, at 1.53 A, has a height of 5.4 in the 
classical RDF, almost 4 times more than the quantum 
result of 1.45. In Fig. 6 we also present the H-H RDF 
obtained by Finney et alP' from neutron diffraction ex- 
periments at 80 K. Note that this curve does not include 
the intramolecular H-H pairs in the original publication. 
The next peak, corresponding to hydrogen pairs in ad- 
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FIG. 6: Hydrogen-hydrogen radial distribution function for 
HDA ice at P — 1 atm. Solid and dashed lines show results 
obtained from PIMD and classical simulations, respectively, 
at T — 75 K. The dashed-dotted line was derived from neutron 
diffraction experiments at 80 K.— 
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FIG. 7: Oxygen-oxygen radial distribution function for HDA 
ice at P = 1 atm. Solid and dashed lines represent results 
derived from PIMD and classical simulations, respectively, at 
T = 75 K. The dashed-dotted line was derived from neutron- 
scattering experiments at 80 K.— 



jacent (H-bonded) molecules appears nearly at the same 
position in the RDF derived from experiment and that 
obtained from PIMD simulations. However, the former 
is higher than the latter. This is the main difference be- 
tween both results, that in general are rather similar. 

Quantum effects are not only associated to hydrogen, 
due to its small mass, but are also observable in the 
oxygen-oxygen RDF of amorphous ice, as displayed in 
Fig. 7. In this figure, the solid line represents g(r) de- 
rived from PIMD simulations at T — 75 K and atmo- 
spheric pressure, and the dashed line corresponds to the 
RDF derived from classical molecular dynamics simula- 
tions at the same conditions of pressure and tempera- 
ture. The classical result is similar to the oxygen-oxygen 
RDF derived in Ref. [U for several interatomic poten- 
tials. For comparison, we also show in Fig. 7 the 0-0 
RDF of HDA ice derived by Bowron et al. from neu- 
tron scattering experiments 4 ^ (dashed-dotted line). As 
in the case of the O-H RDF, we observe a broadening of 
the peaks when nuclear quantum motion is considered, 
as a consequence of the associated atomic delocalization. 
This is particularly observable for the first peak at about 
2.8 A, whose height decreases when quantum effects are 
taken into account. The result derived from PIMD sim- 
ulations is closer to the experimental data than the clas- 
sical one, but the peak derived from the quantum simu- 
lations appears at r = 2.77 A, a distance slightly larger 
than that corresponding to the maximum of the RDF 
derived from experiment (r — 2.72 A). We note, how- 
ever, that the position of the peak derived from another 



neutron diffraction work^i is closer to that obtained in 
our simulations, but its height is somewhat smaller than 
that found here. RDFs of HDA ice have been also derived 
from x-ray diffraction measurements. 7 ? 

It is also interesting to analyze the effect of pressure on 
the shape of the 0-0 RDF of amorphous ice, as it can 
give information on the molecular reorganization in the 
way to amorphous phases with still high density. With 
this purpose, in Fig. 8 we present the 0-0 RDF at 75 K 
for different pressures: atmospheric pressure along with 
P = 1, 3, and 6 GPa. We first observe that the peak 
at about 2.8 A, corresponding to the first coordination 
shell, moves to shorter distances as the pressure is raised, 
in line with a decrease in the 0-0 distance associated to 
the corresponding volume reduction [dV/dP < 0). It is 
also remarkable that the broad feature appearing in the 
RDF around 4 A (at atmospheric pressure) sharpens and 
moves to smaller distances for increasing pressure, indi- 
cating that water molecules in the second coordination 
shell come closer to those in the first shell. This feature 
in the 0-0 RDF coincides with that observed in very- 
high density amorphous (VHDA) icei 14 ' 76 Other features 
of the RDF appearing at larger r also move to shorter 
distances, as observed in the region between 5 and 6 A 
in Fig. 8. As shown above, we do not observe any clear 
discontinuity in the volume as pressure is increased, so 
that according to our results VHDA seems to appear as a 
high-pressure regime of HDA. Although in principle it is 
not evident that the phase obtained by applying pressure 
to HDA ice is the same as the VHDA obtained by iso- 
baric heating, — our results are in agreement with earlier 
simulations, 44 favoring a continuous transition from HDA 
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FIG. 8: Oxygen-oxygen radial distribution function for 
amorphous ice at T — 75 K and several pressures, as derived 
from PIMD simulations. Different lines represent results for 
P = 1 atm, and 1, 3, 6 GPa, as indicated in the labels. 

to VHDA ice. From an experimental point of view, there 
are indications^ pointing in the direction that structural 
differences between HDA and VHDA become less promi- 
nent as pressure increases. In any case, there is some 
evidence against a first order-like nature of the transi- 
tion between HDA and VHDA, suggesting a continuous 
character of the transition^ 



C. Kinetic energy 

In this section we present the kinetic energy of atomic 
nuclei in HDA ice, as derived from our PIMD simulations. 
The kinetic energy, Ek , of atomic nuclei depends on their 
mass and spacial delocalization, so that it can give infor- 
mation on the environment and interatomic interactions 
seen by the considered nuclei. We note that this does 
not occur for classical simulations, since in this case each 
degree of freedom contributes to the kinetic energy by an 
amount that depends only on temperature, i.e., ksT/2. 
A typical quantum effect associated to the atomic motion 
in solids is that the kinetic energy converges at low tem- 
perature to a value related to zero-point motion, contrary 
to the classical result where E/~ vanishes in the limit T — > 
K. Path integral simulations allow us to obtain the ki- 
netic energy of the considered quantum particles. For a 
particle with a certain mass at a given temperature, the 
larger the spread of the quantum paths, the smaller the 
kinetic energy, in line with the expectancy that a larger 
quantum delocalization is associated with a reduction in 
the kinetic energyi 42 i 51 We have calculated here the ki- 
netic energy Ek by using the so-called virial estimator, 
which has an associated statistical uncertainty apprecia- 



FIG. 9: Kinetic energy of hydrogen in ice Ih and amor- 
phous ice as a function of pressure at two temperatures: 75 
K (squares) and 250 K (circles). Open and solid symbols cor- 
respond to ice Ih and amorphous ice, respectively. Error bars 
for the crystalline phase are less than the symbol size. Lines 
are guides to the eye. 



bly lower than the potential energy of the system i 58 i 79 

In Fig. 9 we display Ej. for hydrogen as a function 
of pressure at two temperatures: 75 and 250 K, as de- 
rived from our PIMD simulations of amorphous ice (solid 
symbols). At each temperature, Ek increases slowly as 
pressure rises, corresponding to an overall increase of vi- 
brational frequencies. For comparison, we also present 
results for hydrogen in ice Ih at the same temperatures 
(open symbols). In this case, the kinetic energy increases 
with rising pressure faster than for amorphous ice, in the 
whole region where ice Ih is found to be mechanically 
stable. This is a consequence of the larger quantum de- 
localization of hydrogen in amorphous ice, as compared 
with ice Ih at the same temperature, or equivalently, to 
the presence in the amorphous solid of modes with lower 
vibrational frequency. At atmospheric pressure, Ek for 
HDA ice increases only by about 1.5% from 75 to 250 
K, reflecting the fact that at these temperatures most 
vibrational modes with large hydrogen contribution (i.e., 
libration and stretching modes) are nearly in their ground 
state. Moreover, changing the external pressure modifies 
even less the kinetic energy of hydrogen. In fact, at 75 
K it rises by about 0.15% in the range from atmospheric 
pressure to P = 8 GPa. 

In Fig. 10 we show the pressure dependence of the ki- 
netic energy of oxygen in ice (HDA and Ih) at the same 
temperatures as those presented in Fig. 9 for hydrogen. 
We observe again that the kinetic energy rises with in- 
creasing pressure, and that at a given pressure, Ek for 
oxygen in amorphous ice is smaller than in ice Ih. At 75 
K and atmospheric pressure the kinetic energy of oxygen 
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FIG. 10: Kinetic energy of oxygen in ice Ih and amor- 
phous ice as a function of pressure at two temperatures: 75 
K (squares) and 250 K (circles). Open and solid symbols 
correspond to ice Ih and amorphous ice, respectively. Error 
bars for the crystalline and amorphous solids are less than the 
symbol size. Lines are guides to the eye. 



results to be about 4 times smaller than that of hydrogen, 
as could be expected from the larger mass of the former 
atom. At this temperature, Ek for oxygen rises by about 
22% when a pressure of 8 GPa is applied. This relative 
change is much larger than that found for hydrogen at 
the same conditions (a 0.15%). For increasing temper- 
ature, we find also an important rise in Ek for oxygen, 
e.g., at P = 1 atm we have a change of a 44% from 75 
to 150 K, to be compared with a rise of 1.5% in the case 
of hydrogen. This is indeed a consequence of the much 
larger mass of oxygen, which contributes mainly to vi- 
brational modes with low frequency, and therefore with 
excited states non-negligibly populated at these temper- 
atures. 

We note that the vertical scales in Figs. 9 and 10 are 
different, and what seems to be a very large change in 
the kinetic energy of hydrogen when comparing Ih and 
HDA ices, is in relative terms less than the change in Ek 
for oxygen. In fact, at P = 1 GPa and T = 75 K (close 
to the amorphization pressure of ice Ih), the kinetic en- 
ergy of oxygen in HDA ice is 167 J/mol lower than that 
in ice Ih, to be compared with a decrease of 270 J/mol 
in Ek for hydrogen. In relative terms these differences 
amount to a decrease of 4.5% for oxygen vs. 1.9% for 
hydrogen. This is in agreement with the fact that the 
structural changes associated to ice amorphization affect 
mainly to low-frequency vibrational modes (translational 
modes of the whole water molecule), with a major rela- 
tive contribution of oxygen to the kinetic energy. The 
relative difference between kinetic energies in ice Ih and 
HDA ice decreases as temperature is raised, as observed 



in Figs. 9 and 10 at T = 250 K, mainly because the nu- 
clear motion becomes "more classical," and therefore Ek 
is less sensitive to the environment and actual motion of 
the atomic nuclei under consideration. 

The kinetic energy of hydrogen and oxygen in liq- 
uid water and ice Ih has been studied in detail earlier 
from PIMD simulations, and compared with data de- 
rived from deep inelastic neutron scattering in the case 
of hydrogen42 In that paper, a study of the contribution 
of different vibrational modes to the kinetic energy was 
presented. For HDA ice we find here that, at a given 
temperature, Ek increases as pressure is raised for both 
hydrogen and oxygen. This result is not obvious, since 
the O-H stretching frequencies are known to decrease for 
increasing pressure (their mode Griineisen parameter is 
negativ e 41 ' 80 ). However, the overall contribution of the 
vibrational modes to Ek increases with pressure, since the 
contribution of modes with positive Griineisen parameter 
dominates, as analyzed elsewhere from a quasi- harmonic 
approximation for crystalline ice phases. 80 



D. Bulk modulus 

The compressibility of ice displays peculiar properties 
associated to the hydrogen-bond network. For the crys- 
talline phases of ice, and ice Ih in particular, the com- 
pressibility is smaller than what one could suspect from 
the large cavities present in its structure, which could 
be expected to collapse under pressure before molecules 
could approach close enough to repel each other. This is 
in fact not the case, and for ice Ih the H bonds holding 
the structure are known to be rather stable, as mani- 
fested by the relatively high pressure necessary to break 
down the crystal. 5 Here we present results of PIMD sim- 
ulations for the isothermal bulk modulus of HDA ice at 
different temperatures and pressures, and compare them 
with those derived for ice Ih. 

The isothermal compressibility k of ice, or its in- 
verse the bulk modulus [B = 1/k = -V{dP/dV) T ] 
can be directly derived from PIMD simulations in the 
isothermal-isobaric ensemble. In this ensemble, B can 
be obtained from the mean-square fluctuations of the vol- 



ume, Oy = (V 2 ) — (V') 2 , by using the expressio: 



,81.82 



B 



k B T(V) 



(3) 



kg being Boltzmann's constant. This expression 
has been employed earlier to obtain the bulk mod- 
ulus of different types of solids from path-integral 
simulations i 51 i 82 >^ 

In Fig. 11 we present the bulk modulus of amorphous 
ice as a function of pressure, as derived from our PIMD 
simulations at 75 K (solid circles) and 250 K (solid 
squares). For comparison, we also present results for ice 
Ih at 75 K (open circles). The bulk modulus of amor- 
phous ice is found to increase linearly as a function of 
pressure in the range considered here. From linear fits to 
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FIG. 11: Pressure dependence of the bulk modulus of amor- 
phous ice at T — 75 K (solid circles) and 250 K (solid squares), 
as derived from PIMD simulations. Open circles correspond 
to ice Ih at 75 K. Error bars for the crystalline phase are in the 
order of the symbol size. Dashed lines are guides to the eye. 
The solid line was obtained by numerical differentiation from 
the pressure-density data displayed in the review by Loerting 
and Giovambattista,— adapted from Ref. lofSI . 



the data shown for amorphous ice in Fig. 11, we find a 
slope dB/dP = 7.1(2) and 7.0(2) at 75 and 250 K, re- 
spectively, i.e. both values coincide within the precision 
of our results. This contrasts clearly with the pressure 
dependence of the bulk modulus found for ice Ih (open 
circles in Fig. 11). In fact, for this crystalline phase of 
water one finds that at P > 0.3 GPa, B decreases as 
pressure is raised, until eventually reaching the limit of 
mechanical stability of the phase (spinodal line) and the 
associated amorphization of the material. Thus, at pres- 
sures close to 1 GPa, the bulk modulus of ice Ih is lower 
than that of the amorphous phase, contrary to the re- 
sult obtained at lower pressures. For comparison, we 
also show in Fig. 11 the pressure dependence of the bulk 
modulus of HDA ice at 80 K (solid line) , derived by nu- 
merical differentiation of the pressure-density data given 
in Ref. GJ (adapted fromll). 

At atmospheric pressure and 75 K, we find for amor- 
phous ice B = 9.6(4) GPa, to be compared with a value 
of B = 14.0(3) GPa derived for ice Ih at the same con- 
ditions. This reflects an important increase in the com- 
pressibility of ice upon amorphization, i.e. a lowering of 
the bulk modulus of the material by about a 30%. Note 
that this reduction in the bulk modulus is similar to the 
relative volume change upon amorphization, as discussed 
in Sect. III. 

From our data in the temperature region between 75 
and 250 K, B extrapolates to zero at negative pressures in 
the order of -0.5 to -1 GPa, which gives the limit of me- 



chanical stability of this amorphous material, where the 
solid breaks down giving rise to the gas phase. As dis- 
cussed elsewhere ; 67 ' 84 the possibility of studying a solid 
in metastable conditions, close to a spinodal line is lim- 
ited by the appearance of nucleation events, which cause 
the breakdown of the solid. For atomistic simulations 
such as those employed here, the probability of those nu- 
cleation events at low temperatures is relatively low, and 
the metastable range of the solid that can be explored 
is rather large. In fact, we can go here to conditions 
near the limit of mechanical stability at negative pres- 
sures (limit B — > 0). As T increases the probability of 
nucleation becomes higher, and the accessible pressure 
range in the simulations is reduced. 



V. COMPARISON WITH EARLIER WORK 

Data of earlier simulations of HDA ice have been al- 
ready presented in the previous section along with the 
results of our PIMD simulations. Here we discuss and 
summarize the main similarities and discrepancies be- 
tween our results and those given in some earlier works. 

Simulations of the amorphization of ice Ih using clas- 
sical molecular dynamics were carried out by Tse and 
Klein ; 85 ! 86 who employed the TIP4P intermolecular po- 
tential, and found ice amorphization at pressures around 
1.2-1.3 GPa at temperatures between 80 and 100 K. 
These values are close to the limit for mechanical stabil- 
ity (spinodal pressure, P s ) of ice Ih obtained from PIMD 
simulations, i.e., P s = 1.12 and 1.26 GPa at 50 and 100 
K, respectively^ 7 - 

Seidl et al. 2i carried out detailed classical molecu- 
lar dynamics simulations of HDA ice in the isothermal- 
isobaric ensemble using several force fields. In particular, 
they studied the glass-transition at a pressure of 0.3 GPa, 
and found some indications of a glass-to-liquid transition 
at a temperature around 200 K, which could suggest that 
HDA ice is a proxy of an ultraviscous high-density liquid. 
From our present results, we cannot find any evidence 
of such a transition, and a detailed study of this point 
with PIMD simulations would require at present enor- 
mous computational resources. 

In connection with our work, Tse et alM- studied amor- 
phous ice by classical molecular dynamics simulations 
with the SPC/E potential. They presented 0-0 RDFs 
very similar to those obtained in our classical simula- 
tions (shown in Fig. 7). These authors found a contin- 
uous and smooth increase in the molar volume of HDA 
ice as temperature is raised, without a sharp change in- 
dicating transformation from high-density to low-density 
amorphous phases. This is in line with the results of our 
PIMD simulations shown in Fig. 3. In particular, they 
found at 100 K a molar volume of 15.4 cm 3 /mol, close 
to our classical result (triangles in Fig. 3) and smaller 
than the value derived from our PIMD simulations with 
the q-TIP4P/F potential (v = 16.2 cm 3 /mol). At higher 
temperatures the results by Tse et al. become closer to 
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our PIMD results, and are therefore somewhat higher 
than our classical data. This can be a consequence of 
the differences between the force fields employed in both 
works. 

Earlier studies of path integral simulations of amor- 
phous ice are scarce. Gai et al, 2 - studied structural prop- 
erties of HDA ice at 77 K by path-integral Monte Carlo 
simulations, using the SPC/E potential model, which 
treats the water molecules as rigid bodies. Their simula- 
tions were carried out in the constant volume, canonical 
(NVT) ensemble, and dealt with both H 2 and D 2 ice. 
The O-H RDF obtained by these authors for D2O amor- 
phous ice is similar to that obtained here, and presented 
in Fig. 5. In particular, they found well-defined peaks 
for the first and second coordination shells at about 1.8 
and 3.3 A, respectively. However, for H2O ice they found 
a RDF very different to that found here, in which the 
peaks at 1.8 A (intermolecular O-H bridges) and 3.3 A 
were missing, and had been replaced by a broad feature 
extending from about 2.5 to 4 A. These authors argued 
that the hydrogen bonding network that is present for 
D2O either disappears or is totally mixed with the sec- 
ond nearest-neighbor shell. 

Apart from the constant volume employed by Gai et 
al£^ in their simulations vs the constant pressure em- 
ployed in ours, the main difference between both kinds of 
calculations seems to be the interatomic potential: rigid 
molecules in Ref. [2l| vs flexible molecules in our calcula- 
tions. We do not find, however, a direct explanation why 
the use of rigid molecules should cause a so strong differ- 
ence between the structures of H 2 and D 2 amorphous 
ice, as suggested by the results of Gai et al. 

VI. SUMMARY 

PIMD simulations provide us with a suitable tool to 
analyze effects of nuclear quantum motion in amorphous 
ice at finite temperatures. Here, we have presented re- 
sults of PIMD simulations of HDA ice in the isothermal- 
isobaric ensemble at different pressures and tempera- 
tures. This kind of simulations have allowed us to ob- 
tain structural and thermodynamic properties of this 
metastable material in a large region of pressures, in- 
cluding tensile stresses (P < 0). 

The HDA ice studied here was obtained computa- 
tionally by pressure-induced amorphization of cubic and 
hexagonal ice (Ih and Ic). We observe an important re- 
duction of the volume upon amorphization at a pressure 
of about 1.2 GPa. The resulting ice at atmospheric pres- 
sure and 100 K is 19% denser than its crystalline pre- 
cursor, but it is found to be softer, in the sense that 
its compressibility and thermal expansion coefficient are 
clearly larger. At P = 1 atm and T = 75 K, the com- 
pressibility k of HDA ice is about 50% larger than that 
of ice Ih, and the thermal expansion coefficient a for the 
amorphous solid is found to be 9 x 10~ 4 K _1 , whereas it 



is negative for ice Ih at 75 K. 

We have assessed the importance of quantum effects 
by comparing results obtained from PIMD simulations 
with those obtained from classical simulations. Struc- 
tural variables are found to change when nuclear quan- 
tum motion is considered, especially at low temperatures. 
Thus, the crystal volume, interatomic distances, and ra- 
dial distribution functions suffer appreciable modifica- 
tions in the range of temperature and pressure considered 
here. At 50 K the molar volume of HDA ice is found to 
rise by 0.85 cm 3 /mol (a 6% of the classical value), and 
the intramolecular O-H distance increases by 1.4% due 
to quantum motion. 

The zero-point vibrational motion of atomic nuclei is 
large enough to also change appreciably structural ob- 
servables of the amorphous solid, such as the radial dis- 
tribution function at low temperatures. In fact, from 
PIMD simulations we observe a broadening of the peaks 
in the RDFs, as compared with classical molecular dy- 
namics simulations. For different isotopes we also ob- 
serve a change in the RDFs of HDA ice. In particular, 
the width of the peaks in the O-H RDF is found to de- 
pend on the hydrogen isotope under consideration, but 
the general features of this RDF are basically the same 
for both H 2 and D2O, contrary to earlier path-integral 
Monte Carlo results^ 

At a given temperature, the kinetic energy of both hy- 
drogen and oxygen is found to increase for rising pres- 
sure. This increase is, however, slower than that ob- 
tained in ice Ih. Such an increase is associated to the rise 
in vibrational zero-point energy and an overall increase 
in the vibrational frequencies of the solid. However, the 
intramolecular O-H distance is found to increase as pres- 
sure is raised, with a decrease in the frequency of the 
corresponding stretching modes. In HDA ice the bulk 
modulus is found to increase linearly as a function of 
pressure, in the whole region studied here. At 75 K we 
find dB/dP = 7.1. 

Although quantitative values found by using the q- 
TIP4P/F potential can change by employing other in- 
teratomic potentials, the main conclusions obtained here 
can hardly depend on the potential employed in the simu- 
lations. Quantum simulations similar to those presented 
here can give information on the atom derealization and 
anharmonic effects in other kinds of amorphous ice. An 
extension of this work could consist in studying amor- 
phous ice at still higher pressures, where very-high den- 
sity amorphous ice could be characterized, including nu- 
clear quantum effects. 
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